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In calculating electroweak radiative corrections at two-loop level, one 
encounters Feynman graphs with several different masses on the internal 
propagators and on the external legs, which lead to complicated scalar 
functions. We describe a general analytic-numerical reduction scheme for 
evaluating any two-loop diagrams with general kinematics and general 
renormalizable interactions, whereby ten basic functions form a complete 
set after tensor reduction. We illustrate this scheme by applying it to two- 
and three-point functions. We discuss the treatment of infrared singular- 
ities within this numerical approach. 
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Because of the level of experimental precision attained in measuring the elec- 
troweak parameters at LEP, SLC, and Tevatron, a full two-loop analysis is desirable. 
This will be even more necessary in view of the precision envisioned at future colliders 
such as the LHC, NLC, and the GigaZ. 

While several two-loop quantities are already included in standard electroweak 
fitting programs such as ZFITTER they are often obtained within certain ap- 
proximations where one can neglect certain masses or can perform a mass expansion 
of Feynman graphs. These approximation techniques, while ingenious, were used 
because the exact Feynman integrals typically lead to complicated scalar functions 
which often cannont be evaluated analytically, in a closed form, in terms of usual 
special functions. During the past decade, the existing work on massless and mas- 
sive two-loop graphs made it clear that for the general mass case, a certain 
amount of numerical work is unavoidable. 

Here we discuss the status of a hybrid, analytical-numerical approach to two- 
loop radiative corrections with arbitrary masses. The aim is to treat any two-loop 
graph, of any topology, by using the same algorithm, so that the general recipe can 
be encoded in a computer program. Such an approach was developed in ref. [Q, and 
was successfully applied to several physical processes 

At the center of this approach is the introduction of a set of ten basic functions, 
hi — hio, defined by the following one- dimensional integral representations: 
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Here we use the following notations: 
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The evaluation of these ten functions is best done by numerical integration. 
By using an adaptative deterministic numerical integration algorithm, these one- 
dimensional integrals can be calculated fast and very precisely. By numerical inte- 
gration, we plot these ten basic functions in figure 1 for a range of their kinematic 
variables. 

Within the method we discuss here, any two-loop diagram with arbitrary masses is 
first reduced to multi-dimensional scalar integrals involving these ten basic functions 
hi. In ref. [§] we have shown that this can always be done for any two-loop diagrams 
occuring in renormalizable theories. We note in passing that for the case of non- 
renormalizable theories this set of ten functions hi in general needs to be extendeed 
to include additional functions. 
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Figure 1: Plots of the ten basic functions /ij(m^, m^, m^; —k'^) as a function of the external 
momentum variable k"^. The plots given here are for = = m| = 1. 





Figure 2: Expressing generic massive two-loop Feynman diagrams as integrals over sunset- 
type functions. 

This reduction procedure starts by introducing Feynman parameters in the orig- 
inal Feynman graph in order to relate the integral to sunset-type integrals. This is 
illustrated in figure 2. 

At his point, the Feynman graph is written as a multi-dimensional integral over 
tensor integrals of the following type: 



In the following step, these tensor integrals are decomposed into scalar integrals, 
whereby the Lorentz structure is constructed from the vector k'^ and the metric tensor 
g'^'^. This can be carried out systematically by decomposing the loop momenta p an 
q into components parallel and orthogonal to k^. For instance, the decomposition of 
all tensor integrals up to rank three is the following: 
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In the formulae above, a loop integration / (T'p d"q is understood, and we used the 
following notations: 

[211] = [ip + ky + mlf{q' + ml){r' + ml) , r'^'^ = g^^ - (7) 

In ref. [^ we have shown that all the scalar coefficients iAj involved in this tensor 
decomposition are directly expressible in terms of the ten basic functions hi, up to 
trivial one-loop tadpole integrals. In all the formulae above, we have considered only 
integrals with a special combination of propagator powers, namely [211] of eq. 6. 
Where integrals with higher powers are needed, they can be obtained directly by 
mass diferentiation. The only two- loop combination with lower power, [111], can be 
obtained from [211] by a recursion formula obtained by partial integration |^. 

After performing these steps, the Feynman graph is decomposed into scalar inte- 
grals expressed essentially as multiple integrals over hi functions, plus trivial one-loop 
tadpole-type contributions. All necessary formulae to perform this reduction are given 
in ref. [^. They were encoded into computer algebra programs for automatizing the 
reduction. 

Once this standard integral representation is obtained for all Feynman graphs in- 
volved in a physical process, the final step consists in a numerical multi-dimensional 
integration of these expressions. The numerical integration uses an adaptative de- 
terministic algorithm, similar to the numerical integration for the hi functions. This 
ensures an efficient and precise evaluation of the integrals. 

Several physical calculations have been performed so far by using this method. 
As a two-point example to test the reduction algorithm and the reliability of the 
numerical integration, in figure 3 we show the mixed electroweak-QCD Feynman 
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Figure 3: The two-loop Feynman graphs which contribute to the b-mass dependent correction 
of 0(0^5^) to the top self-energy. Only the counterterm diagrams are shown which are 
needed for subtracting the infinities of the imaginary part of the self-energy, which gives the 
0{as) correction to the t — > TV + 6 decay. 
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160 
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1.127 
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1.402 


1.553 


1.712 




-.092 


-.104 


-.117 


-.132 


-.149 



Table 1: The 0{as) correction to the top decay t ^ W -\-b as obtained from the imaginary 
part of the two-loop top self-energy of figure 1, integrated numericaUy. We took Gp = 
1.16637 • 10"^ GeV-2, mw = 80.41 GeV, = 4.7 GeV, and as{mt) = .108. 

graphs which contribute to the top quark self-energy at two-loop. By calculating the 
self-energy function S(p-7) = Si(p- 7) -I-75 ■ S^5(p- 7) at two- loop, from its imaginary 
part one can extract the top decay width up to 0(as), as Ft = 2 ■ JmSi(p ■ 7 = mt). 
Since this correction is known in an analytic form, this provides a good check on our 
two-loop algorithm. The results for the correction to the width, obtained from the 
imaginary part of the two- loop self-energy, are given in table 1. They agree with the 
existing analytic results. 

As a three-point example, in figure 4 we show the diagrams which contribute to 
the top-dependent decay process Z —>■ bb. 

A point which deserves special attention in this case is the presence of IR di- 
vergences. Because the general formulae of the hi functions are derived for general, 
finite masses, IR divergences in our method are not automatically extracted as poles 
in the space-time regulator, as is costumary in QCD calculations. IR divergences 
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Figure 4: Two-loop three-point diagrams contributing to Z — > 66 at 0{asg^). 

in our approach usually appear as end-point singularities in the Feynman parameter 
integration over hi functions, and therefore require a special treatment. 

For both cases at hand - the top quark self-energy and the Z —>■ bb decay - one 
possible approach is to use a mass regulator for the gluon. While in general this is 
not possible for non-abelian theories becaus it does not preserve the Slavnov- Taylor 
identities, in the particular case of these mixed electroweak-QCD corrections this is a 
correct procedure because at this order the IR structure is the same as in the abelian 
case. 

Another approach is to extract the IR structure of the graphs analytically before 
numerical integration. This can be done in the form of one-loop integrals which can 
be handled separately in an analytical way, by usual dimensional regularization. This 
is illustrated in fig. 5. Once the IR divergences are extracted in the form of one- loop 
integrals, the two-loop integration can be carried out numerically. 

We give in table 2 numerical results for all two-loop Feynman graphs involved in 
this process. The numerical results are after the extraction of the UV poles. The IR 
singularities are subtracted as shown in figure 5. The numerical integration accuracy 
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Figure 5: Extracting the infrared divergent pieces of the two-loop diagrams analytically. The 
infrared divergency of the two-loop diagram is the same as the infrared divergency of the 
product of the two one-loop diagrams obtained by "freezing" the common line in the loop 
momenta integration. 
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Table 2: Numerical values for the two- loop diagrams shown in figure^ Vi-Vio are the sums 
of the corresponding W and (j) exchange graphs. An overall color and coupling constant 
factor of 27^(1—75)0^(^3/12 cos 6w) is understood. The UV and IR divergences are removed 
as discussed in the text. 



is 10-3. The evaluation of a total of 78 Feynman graph evaluations with this precision 
required 100 hours computing time on a 600 MHz Pentium machine. 

To conclude, we developed an algorithm for the tensor reduction of massive two- 
loop diagrams. It applies in principle to any massive two-loop graph, and it can be 
automatized in the form of a computer algebra program. The tensor decomposition 
algorithm results in a multi-dimensional integral over a set of ten basic functions hi, 
which are defined in terms of one- dimensional integral representations. We described 
the numerical methods which we used for carrying out the remaining integrations. 

We have shown how these techniques work in the case of two realistic calculations 
of mixed electroweak-QCD radiative corrections. The first example is the two-loop 
top quark self-energy from which the 0{as) correction to the top quark decay width 
can be extracted and compared with the analytical result. The second example is 
a three-point calculation involving all two-loop diagrams which contribute to the 
top-dependent decay process Z bb. Thus we have shown that the techniques 
we described can be used in realistic calculations, where several internal mass and 
external momenta scales are involved. 

This approach works for any such combination of kinematic variables, apart from 
possible infrared complications. In the context of the Z ^ bb example, we discussed 
the analytical separation of the infrared divergencies. Within our two-loop methods, 
if a process involves infrared singularities, these have to be dealt with in a special 
way because the numerical nature of our methods. 
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